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Abstract 

A new splitting scheme is proposed for the numerical solution of the time- 
dependent, incompressible Navier-Stokes equations by spectral methods. A 
staggered grid is used for the pressure, improved intermediate boundary 
conditions are employed in the split step for the velocity, and spectral 
multigrid techniques are used for the solution of the implicit equations. 


Research for the second author was supported by the National Aeronautics 
and Space Administration under NASA Contract No. NAS1-17070 while he was in 
residence at ICASE, NASA Langley Research Center, Hampton, VA 23665. 


i 




I. INTRODUCTION 


Numerical simulations of transition and turbulence are playing an 
increasingly important role m the investigation of the basic physics of these 
strongly nonlinear, three-dimensional phenomena. The most efficient direct 
simulation codes are those which use spectral methods in simple geometries 
such as the periodic box ([1], [2]), the channel ([3] - [7]), the pipe ([8], 
[9]), the parallel boundary layer ([7], [8], [10]), and the region between two 
infinite, concentric cylinders [11]. All of these algorithms require periodic 
boundary conditions in at least two directions and thus are unsuitable for 
such more realistic problems as the true, non-parallel boundary layer and the 
flow between two finite, concentric cylinders. Moreover, in most applications 
the viscosity must be treated implicitly, at least in the non-periodic 
directions, and these algorithms are quite limited in the nature of the 
viscosity that they can handle. The most flexible of these algorithms [7] can 
treat implicitly a viscosity which varies in time and in the one non-periodic 
direction. 

The limitation to at most one non-periodic direction arises from the 
solution schemes used for the implicit systems of equations in the algorithms, 
typically a Poisson equation for the pressure and Helmholtz equations for the 
velocity components. These are usually solved by direct methods which first 
diagonalize the equations in the periodic directions and then resort to 
special efficient solution procedures which apply only when there is at most 
one non-periodic direction. These spectral methods invariably use Fourier 
series in the periodic directions and Chebyshev or some other Jacobi 
polynomials in the non-periodic directions. The Fast Fourier Transform (FFT) 
can be used to diagonalize the equations in the periodic directions provided 
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that the viscosity (or at least the part of it which is treated implicitly) is 
independent of the periodic coordinates. In a limited number of cases an 
efficient solution procedure exists for the remaining, non-periodic direction. 
Thus, in addition to the obvious need for new, more general spectral 
algorithms for problems with more than one non-periodic coordinate, there is 
also a need for new algorithms even for fully periodic problems with strongly 
position-dependent viscosity. 

The way to achieve greater generality is to abandon direct solution 
techniques in favor of iterative methods. The algorithm described in [7] does 
use iterative techniques, but only in one of the directions. This method does 
not use splitting and the implicit equations are solved together as a system. 
A preconditioned iterative scheme is at the core of this method. The 
effective iteration matrix has complex eigenvalues. The extension of this 
method to the more general cases of interest here is not straightforward. The 
use of splitting would lead to implicit equations which are positive definite 
and much more straightforward to solve by iterative methods. However, Marcus 
[11] has clearly pointed out serious errors that can arise in splitting 
methods. In this paper we describe a new type of splitting in which these 
errors are absent. Moreover, we discuss efficient iterative solution of the 
implicit systems by spectral multigrid techniques ([12] - [16]). 
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II. SPLITTING SCHEME 

The incompressible Navier-Stokes equations can be written 
_u = u x a) + V*(pVu) - VP in 

in n 
on 3ft 


V«u = 0 


u = g 


( 1 ) 


u_(x,0) = Uq(£) in ft, 

where _u_ = (u,v,w) denotes the velocity, a = V x u the vorticity, 

P = p + -i- |u| the total pressure, p the viscosity, ft the interior of the 
domain, and 3ft its boundary. In some problems the above system should be 
supplemented with an equation for the temperature, and the viscosity may 
depend upon position x and time t. An empirical formula for the viscosity 
is given in [17] (and repeated in [7]) in terms of the temperature. A 
complete description of the equations governing the parallel boundary layer is 
given in [7] . 

The solution algorithm is based upon the following splitting scheme for 
advancing from t = t n to t = t n+ ^ : 

* * * * 

£ = ii x (o + V»(pVii ) in ft 

u = g on 3ft (2) 

•k 

U ( X , t ) = U< X , t ) 


in ft, 
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** ** 

_u fc = - V P in fi 

** 

V*_u = 0 in n 

ii'k A A 

u_ • n = _g • n on (3) 

** <* ** A 
t£ t • t = - V P • t on 3!1 

** * * * 
u (x,t ) = u (x,t ) in Q, 

where Eq. (2) is integrated from t n to the intermediate time t ; then Eq. 

•J(f 

(3) is integrated from t to t n+l and ii(x., t^ + ^) * s identified with 

** a A 

ix (xc,t n+ j). Here n denotes the unit normal to the boundary and t a unit 

tangent vector at the boundary. Note that in the velocity step all of the 

velocity components are specified at the boundary, whereas in the pressure 

step only the normal velocity component is specified on the boundary, and the 

differential equation itself is used for the tangential velocity components 

there. The boundary conditions for each step are determined by the condition 

that the individual step be well-posed. At the end of a full step the flow 

will be divergence-free, but there will be a slip velocity at the boundary. 

"k n 

If the boundary conditions on the intermediate, velocity step are u = g or 
u* = g n+ ^ , then the slip velocity at the end of the full step will be 
0(At). This can be reduced by resorting to intermediate boundary conditions 
of the type suggested in [18]. We set, 
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U* • J = (g n + At(g£ + VP n ) ) . T 


* * 
u • n 


(/ + At gj) • n. 


(4) 


o 

This reduces the final slip velocity to 0(At ). Higher-order corrections 
are also possible. 

There are two key elements of this splitting scheme which distinguish it 
from the standard splitting scheme for spectral algorithms [3], whose problems 
have been documented by Marcus [11]. The first is the use of a staggered grid 
for the pressure [7], and the second is the use of intermediate boundary 
conditions, i.e., _g_ , which are not the desired boundary conditions on ju_, 
i.e., _g_, but rather are derived from the principles developed by LeVeque and 
Oliger [18] and applied by Kim and Moin [19] to a two-dimensional finite 
difference algorithm. 

The advantages of the staggered grid for finite difference and finite 
element schemes are well-known. Principal among them is that the pressure 
step does not require any numerical boundary conditions on the pressure. The 
usual procedure on a non-staggered grid is to reduce Eq. (3) to a Poisson 
equation for the pressure and to use the normal momentum equation to obtain a 
Neumann boundary condition on P. However, the incompressible Navier-Stokes 
equations with boundary conditions solely on the velocity form a well-posed 
mathematical problem. The imposition of pressure boundary conditions is an 
overspecification of the problem and can lead to numerical instability or 
inconsistent results. Such difficulties are far more likely to be observed in 
solutions obtained via spectral methods rather than via finite difference 
methods. The former are notoriously sensitive to boundary conditions — see 
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Figure 2 of [20] for a dramatic example of instability arising in a spectral 
calculation with overspecified boundary conditions — and produce solutions that 
are sufficiently accurate so that even minor inconsistencies are noticeable — 
see [11] for an example of a spectral calculation which gives an incorrect 
result because of the pressure boundary conditions. Kleiser and Schumann [4] 
and Marcus [11] have devised a method for imposing a consistent pressure 
boundary condition on the Poisson equation, but this applies only to 
incompressible flow with constant viscosity, and not to the more general 
problems of interest here. 


III. DISCRETIZATION 

For simplicity we will describe the spectral discretization for a two- 
dimensional channel flow problem with periodic boundary conditions in x and 
no-slip boundary conditions in y on £2 = [0,2ir] x (-1,1). The spectral 
method is based upon Fourier series in x and Chebyshev polynomials in y. 
The collocation points in these directions are 


x. = 


2iri 


i N 


i = 0,l,...,N x -l 


y. = cos 
J N 

y 


j = 0 , 1 , • • • , N 


(5) 


y v = cos J - U+Vi 

y 3+ V 2 N y 


j = 0, 1 , • • • >N y -l ■ 


The velocity is defined at (x^,yj) and the pressure at (x^,yj + y ) . The 
momentum equation is enforced at the interior velocity nodes (cell edges) and 
the continuity equation at the pressure nodes (cell centers). 
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The velocity step is discretized in the usual spectral collocation manner. 
A Fourier transform in x simplifies the pressure step. Define the Fourier 
components by 


"k ( V -TT X) 

X £=0 


-ikx 

"(Xj.yjJe 1 


( 6 ) 


for k = -N x /2, -N x /2+l,»*», N x /2-l. Let Cq and C + be the matrices which 
represent the operations of computing Chebyshev coefficients from function 
values at the velocity nodes and the pressure nodes, respectively; let D 
represent the differentiation operator in terms of the Chebyshev coefficients, 
and let iK represent the differentiation operator in terms of the Fourier 
coefficients. (K is a real diagonal matrix and D is a real, strictly upper 
triangular matrix.) 

The discrete divergence operator is 


Vu_ = (C; 1 C Q ) [ iKu + Cq L D C q v], (7) 

(The operator c” 1 Cq represents interpolation from the cell edges to the 
cell centers.) This may be decomposed into a part which depends on the normal 
velocity component at the boundary, denoted by 8_u and the remaining part, 

a 

denoted by ( V - B)n. The discrete gradient operator is 

GP = (C" 1 C + )(iKP, C^ 1 DC + P). 


A A 

Letting Q = AtP, the fully discrete version of Eq. (3) is 


( 8 ) 
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'n+1 


u 


A J» A 

u - 6Q 


at interior cell edges 


"n+1 


u 


u* - (G Q) 


( 9 ) 


'n+1 


^n+K 
(g > 


at y = ±1 


p‘ u n+1 - 0 


at cell centers, 


( 10 ) 


These combine to yield 


flu n+1 = ( V - B)(u* -GQ) + 8g n+1 = 0 


(ID 


or 


IQ = F 


( 12 ) 


with 

LQ = -K 2 Q + C" 1 DC ± C" 1 DC + Q 


(13) 


F 


Vu 


B(g 


n+1 



(14) 


where is Cq with its first and last rows replaced with zeroes. The 
pressure is computed from Eq. (12) and the velocities are then adjusted via 
Eq. (9). The equation for the pressure is diagonal in the wavenumbers k. It 
is singular for k = 0, but this merely reflects the indeterminacy of the 
pressure to within an additive constant. Note that no pressure boundary 
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condition is included in Eq. (12). The right-hand side, however, contains the 
desired boundary conditions on the normal velocity. 

Traditional time discretizations for the Navier-Stokes equations use 
explicit Adams-Bashforth or Runge-Kutta methods for the explicit terms and 
Crank-Nicolson for the implicit terms. We use the low-storage, third-order 
Runge-Kutta/Crank-Nicolson scheme detailed in [7] for the velocity step, with 
a backward Euler pressure correction applied after each stage of the Runge- 
Kutta method. 

The velocity step contains Helmholtz equations for the velocity 
components. If the viscous term in the periodic directions is explicit, this 
equation will have the form (in Fourier space) 

( ~ 2At) i + !y (^ ( y» t) i r) = V (15) 


if it is semi-implicit, 


* 2 — a a /— 

(-2*0^ - k w(y,t)u k + |y I u(y,t) w - 


= f, 


(16) 


and if it is fully implicit, it will have the form (in physical space) 


. 3 


3u ^ 


(- 2A t)u + ^ (y(x,y,t) ^ ^ (ji(x,y,t) = f , 


ay 


(17) 


N -1 

X 


u(y 


>t) = Xw »*Cx 1 ,y, t ). 

x , „ 


i=0 


where 
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Further details on spectral collocation for a similar, but unsplit method 
are available in [7]. The z direction is treated in a manner analogous to 
the x or y directions here, depending upon whether the boundary conditions 
in the third direction are periodic or non-periodic. Of course, if the y 
direction is periodic as well, then it is treated the same as the x 
direction. Spectral collocation methods are readily applied in conjunction 
with mappings. In particular, for the parallel boundary layer problem the y 
direction is semi-infinite. Both algebraic [7] and exponential [10] mappings 
are feasible. 


IV. MULTIGRID ASPECTS 

The role of spectral multigrid (SMG) methods in this spectral collocation 
algorithm for the incompressible Navier-Stokes equations is to accelerate the 
iterative solution of Eq. (12) and either Eq. (15), Eq. (16), or Eq. (17) at 
each time step. We refer the reader to [12] and [13] for the fundamentals of 
SMG and to [21] for a general background on multigrid methods. The trigono- 
metric interpolation techniques described in [13] for the prolongation and 
restriction operators apply directly to the current problems. As usual, the 
relaxation procedures must be tailored to the specific problem. 

In some fully periodic problems, the relaxation may be explicit, but in 
most cases it is essential to precondition the relaxation with a finite 
difference approximation to the differential equation [22]. If the 
preconditioning is used in only one coordinate direction, then multigrid 
methods are not really necessary, for the finite difference approximation may 
be inverted quickly and relaxation on a single grid has a spectral radius 
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independent of the size of the problem. However, if preconditioning is used 
in more than one direction, multigrid methods are important because they allow 
incomplete solutions of the finite difference operator to serve as effective 
and efficient preconditioners. 

Since the FFT may be used for both Fourier series and Chebyshev polynomial 

manipulations, the asymptotic operation count for evaluating the explicit 

terms in both the velocity and pressure steps is 0(N x N^ logCN^N^)) in two 

dimensions. A single iteration on the implicit systems costs 

0(N x Ny log(N x Ny)) operations for the residual evaluation plus whatever work is 

required to invert the finite difference preconditioning operator. A complete 

inversion of the finite difference operator has an operation count which is 

generally of larger order than that for the residual evaluation, whereas an 

incomplete inversion costs 0(N„N„). The most desirable situation would be to 

a y 

have the number of iterations required for an acceptable solution to the 

implicit equations be independent of N„ and N while using incomplete 

a y 

inversions of the finite difference operator. This is not achievable for 
single-grid iterative schemes ([13], [15]). With SMG methods, however, the 
ideal situation is nearly achievable, with the number of iterations growing as 
some small fractional power of N„N„. 

a y 

SMG methods have yet to be applied to fully three-dimensional Poisson 
problems. We anticipate that as in two dimensions, successful smoothing and 
grid coarsening strategies for three-dimensional SMG methods will be the 
appropriate modifications of strategies for three-dimensional finite 

difference multigrid methods (FDMG). On the more difficult three-dimensional 
problems FDMG methods have had to resort to alternating plane relaxation (with 
coarsening in all three directions) or to alternating line relaxation (with 
coarsening in only two directions). 
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Fully Periodic Problems 

For this class of problems the pressure equation is quite simple, and it 
can be solved easily by a direct method since it is a diagonal system in 
Fourier space. The Fourier SMG techniques developed in [12] and refined in 
[16] will solve even Eq. (17) very efficiently. As long as the grid size is 
the same in all directions, Richardson relaxation combined with residual 
averaging and position-dependent relaxation parameters will achieve an order 
of magnitude reduction in the residual per fine grid relaxation [16]. On more 
general grids, SMG is best applied in conjunction with incomplete finite 
difference preconditioning. Brandt, et al. [16], have demonstrated that 
alternating direction Zebra relaxation using the finite difference 
preconditioner is effective on a vector processor. 

Problems with One Non-Periodic Direction 

In this case as well, the pressure equation can be solved efficiently on a 
single grid, since a complete inversion of the preconditioning operator 
requires just the solution of independent tridiagonal systems. Simple 
estimates of the type used by Orszag [22] suggest that the condition numbers 
of the preconditioned systems, i.e., the ratios of the largest to the smallest 
eigenvalues of the effective iteration matrices, should be a decreasing 
function of the Fourier wavenumber, with a maximum of 2.5 for k = 0 and a 
minimum of 1 as k+ ®, A residual reduction by at least a factor of 3 is 
easy to achieve, e.g., by a Chebyshev semi-iterative scheme [22]. The minimum 
residual scheme used in [7] is nearly as effective and requires no relaxation 


parameters . 
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For the velocity step, only the fully implicit viscous time discretization 
requires SMG. Both Eqs. (15) and (16) can be solved efficiently by the same 
techniques used for the pressure equation. For Eq. (17), however, SMG is 
practically a necessity, especially for large problems. Preconditioners for 
Fourier-Chebyshev SMG methods are already available, e.g., ADI relaxation [13] 
as well as incomplete LU decomposition and various Zebra schemes (T. Phillips, 
private communication). 

Problems with Two Non-Periodic Directions 

Here SMG methods are advisable for all the implicit equations because of 
the rapid increase in expense of direct methods for inverting the 
preconditioning operator. Incomplete LU decomposition combined with non- 
stationary Richardson relaxation has already proven successful for problems 
with two Chebyshev directions [13] as have ADI methods [14], Although the 
smoothing rates of the ADI schemes are not as fast as those of the incomplete 
LU schemes , the ADI methods may nonetheless be superior on current 
supercomputers because they vectorize better and require less auxiliary 
storage. 


V. NUMERICAL EXPERIENCE 

The novel aspects of the proposed spectral algorithm are the staggered 
grid in a splitting scheme, the intermediate boundary conditions given in Eq. 
(4), and the SMG solution of the pressure equation. The remaining components 
of this algorithm have been tested extensively elsewhere. 
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The results of the unsplit spectral collocation algorithm described in [7] 
will be used to check the soundness of the present discretization. The first 

test problem will be one of the examples used in [7] to verify the unsplit 

algorithm for channel flow at a Reynolds number of 5000, namely the three- 
dimensional Tollmein-Schlichting wall mode listed there in Table I. Periodic 
boundary conditions are applied in x and z and no-slip boundary conditions 
in y. The grid is 4 x 48 x 4 and the time-step At = 0.10. The initial 

condition consists of the mean flow plus a very small multiple of the linear 

eigenfunction — it disturbs the mean flow by at most 0.001% in the x 
direction. The criterion is how well the computed growth rate matches the 
linear growth rate at t = 50. 

The results appear in Table I. The linear growth rate is -0.076227. Note 

the extremely high accuracy of the unsplit code — even after more than two full 

periods of the Tollmein-Schlichting wave, there is more than five-digit 

accuracy. Results are given in the table for the split code under two 

different boundary conditions. The zero-order boundary conditions are that 
* 

g_ = 0, and the first-order boundary conditions are given by Eq. (4). The 
first-order boundary conditions produce much better results than the zero- 
order boundary conditions. Indeed, the unsplit code performs only slightly 
better than the split code with first-order boundary conditions. 

A second comparison of the split and unsplit codes was performed for a 
Reynolds number of 1500 on a 32 x 64 x 32 grid. The initial conditions were 
essentially the same as those used in reference 23. The time history of the 
energy in several of the horizontal Fourier harmonics is provided in Figure 1 . 
The shear component 9u/9y is given in Figure 2 at t = 15 for the plane 
z = 0. The two calculations agree to better than 0.1%. 
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It might appear preferable to use the same intermediate boundary 
conditions on the normal velocity component as are used for the tangential 
velocity. However, calculations performed with such first-order boundary 
conditions prove to be unstable. This channel flow problem is yet another 
example of the great sensitivity of spectral methods to the boundary 
conditions. Calculations performed with a similar split code which uses 
finite differences rather than Chebyshev polynomials in the normal direction 
display no problems with first-order boundary conditions on all three-velocity 
components. Indeed, calculations performed with the split Fourier-finite 
difference code exhibit a substantial improvement when the full first-order 
intermediate velocity boundary conditions are used in place of the usual zero- 
order ones: for a calculation similar to the one shown in Figure 12 of [7], 
but on a 16 x 100 x 16 grid, the split code agrees with the corresponding 
unsplit code at t = 30 to five significant digits with the first-order 
intermediate boundary conditions, whereas only three-digit agreement is 
achieved with the zero-order ones. 

Several alternative intermediate boundary conditions on the normal 
velocity are possible, and these are currently under investigation to 
determine their accuracy and stability for the fully spectral algorithm. 

As noted above, the pressure equation for this problem may be solved 
acceptably fast with a single-grid relaxation scheme such as the minimum 
residual method. Some preliminary calculations have been made with a SMG 
scheme employing minimum residual relaxation. They do indicate that the fine 
grid relaxations are more effective when combined with coarse grid 
relaxations, e.g., upon returning to the fine grid the residual is initially 
reduced at over twice the rate that it would have been had the calculation 
simply remained on the fine grid. 
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VI. CONCLUDING REMARKS 

The techniques discussed here allow spectral methods to be applied 

reliably and efficiently to a wider range of time-dependent Navier-Stokes 

calculations than heretofore possible. They allow solutions of fully periodic 

problems of greater generality than even those covered by the algorithm 

presented in [2]. They appear to be just as efficient as the best available 

algorithms for parallel channel flow. Together with [7] this is the only 

spectral algorithm suitable for the heated, parallel boundary layer, as well 

as the only spectral algorithm for the usual parallel boundary layer and the 

infinite concentric cylinder problem with a asymptotic operation count as low 

as of N N N log(N N N )}. Moreover, the iterative schemes employed in this 
v x y z ° x y z ' 

algorithm are more robust than those in [7]. 

This algorithm is directly applicable to the finite concentric cylinder 
problem. The collocation points will be automatically clustered near the 
boundaries, where the most resolution is required. It can also be applied to 
the non-parallel boundary layer. Some modifications are desirable because the 
usual collocation points are not well-suited to the resolution requirements in 
the streamwise direction. This difficulty can be alleviated to some extent by 
a streamwise coordinate mapping. A more desirable approach appears to be a 
combination of the best aspects of this method with those of the spectral 
element method ([24], [25]), which in its current form still relies on direct 
solution techniques for the implicit equations. 
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Effect 

Table I. 

of Intermediate Boundary 

Conditions 

Code 

u i ^ calc. 

“J error 

Unsplit 

-.076234 

-.000007 

Zeroth-order 

-.075849 

.000378 

Split 



First-order 

-.076244 

-.000017 

Split 







Figure 1. Harmonic histories for two channel flow simulations. The harmonics 
are labelled by their respective streamwise (x) and spanwise (z) 
wavenumbers; (a) unsplit code; (b) split code. 




X 


Figure 2. 


Vertical shear 3u/9y at t = 15 in the plane z = 0; 
(a) unsplit code; (b) split code. 
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